#include "includesT.h"
#include "pocsT.h"
#include "ScanDataT.h"

using namespace std;

void prepare(ScanData & s)
// function to be executed after the data has been entered in the DataScan object. After this function, the data is ready to start iterations
{
	int i,j,k;

	// Set line mask
	for(k=0; k<s.Nz;k++){
		for (j=0; j<s.Ny; j++){
			for(i=0; i<s.Nx;i++){
				if(s.DBPmask[i][j][k] && s.APmask[i][j][k]) s.Lmask[j][k] = true;
				//if(s.DBPmask[i][j][k]) s.Lmask[j][k] = true;

			}
		}
	}
	
	// Set forward project
	P4prepare(s);	

	// Perform backprojection to set known Hilbert (DBPdata)
	backproj_ibp(s);

	// zero mean test
	/*Float mean=0;
	int N=0;
	for(k=0; k<s.Nz;k++){
		for (j=0; j<s.Ny; j++){
			for(i=0; i<s.Nx;i++){
				if(s.DBPmask[i][j][k]) {mean += s.DBPdata[i][j][k]; N++;}
			}
			mean = mean/N;
			N=0;
			for(i=0; i<s.Nx;i++){
				if(s.DBPmask[i][j][k]) s.DBPdata[i][j][k] = s.DBPdata[i][j][k] - mean; 
			}
			mean=0;
		}
	}*/

}

void P4prepare(ScanData & s){
// This function assumes that the center of the pi-line is at x=0.
	int i,j,k,pi, pf;
	Float x=0,y,z,U,u,v,bterma1, bterma2, btermb1, btermb2, ta, tb,ta1,ta2,tb1,tb2,h;
	for (j=0; j<s.Ny; j++){
		y = s.ymin + s.dy * j;
		ta = ShiftPosition( asin(y/s.R) );
		tb = ShiftPosition( PI - ta );
		pi = (int) PositionIndex(s ,ta,true);
		pf = (int) PositionIndex(s ,tb,false);
		for(k=0; k<s.Nz;k++){
		//for(k=255;k<=256;k++){ // only for central slice
			//Sum the a priori data along the pi-line
			Float APsum = 0;
			for(i=0;i<s.Nx;i++){
				if(s.APmask[i][j][k]) {
					APsum += s.dx*s.APdata[i][j][k];
				}
			}


			z = s.zmin + k*s.dz;
			
			U = s.R - y*s.sint[pi]; 
			u = s.D/U*(-y*s.cost[pi]);
			v = s.D/U*z;
			bterma1 = DFint(s, s.data, pi, u, v)*0.5; //*s.D/sqrt(s.D*s.D+u*u+v*v);

			U = s.R - y*s.sint[pi-1];
			u = s.D/U*(-y*s.cost[pi-1]);
			v = s.D/U*z;
			bterma2 = DFint(s, s.data, pi-1, u, v)*0.5; //*s.D/sqrt(s.D*s.D+u*u+v*v);

			U = s.R - y*s.sint[pf]; 
			u = s.D/U*(-y*s.cost[pf]);
			v = s.D/U*z;
			btermb1 = DFint(s, s.data, pf, u, v)*0.5; //*s.D/sqrt(s.D*s.D+u*u+v*v);

			U = s.R - y*s.sint[pf+1]; 
			u = s.D/U*(-y*s.cost[pf+1]);
			v = s.D/U*z;
			btermb2 = DFint(s, s.data, pf+1, u, v)*0.5; //*s.D/sqrt(s.D*s.D+u*u+v*v); 

			ta1 = IndexPosition(s,pi);
			ta2 = IndexPosition(s,pi-1);
			tb1 = IndexPosition(s,pf);
			tb2 = IndexPosition(s,pf+1);

			s.Fdata[j][k] = bterma1*(ta-ta2)/(ta1 - ta2) + bterma2*(ta-ta1)/(ta2 - ta1) + btermb1*(tb-tb2)/(tb1 - tb2) + btermb2*(tb-tb1)/(tb2 - tb1) - APsum;
		}
	}
}


Float DFint(const ScanData & s, const vector< vector< vector<Float> > > & Ddata , int p, Float u, Float v){
	// Performs bilinear interpolation
	Float wu1,wu2, u1,u2, wv1, wv2, v1,v2;
	int n1,n2,q1,q2;
	
	n1 = (int) floor((u-s.umin)/s.du);
	n2 = (int) ceil((u-s.umin)/s.du);
	u1 = s.umin + s.du*n1;
	u2 = s.umin + s.du*n2;
	if(n1==n2){
		wu1=0.5;
		wu2=0.5;
	}else{
		wu1 = (u - u2)/(u1-u2);
		wu2 = (u - u1)/(u2-u1);
	}

	q1 = (int) floor((v-s.vmin)/s.dv);
	q2 = (int) ceil((v-s.vmin)/s.dv);
	v1 = s.vmin + s.dv*q1;
	v2 = s.vmin + s.dv*q2;
	if(q1==q2){ // point in middle
		wv1=0.5;
		wv2=0.5;
	}else{ // linear interpolation weights
		wv1 = (v - v2)/(v1-v2);
		wv2 = (v - v1)/(v2-v1);
	}
		
	if(q1<0 || q1>=s.Nv){
		if(n1<0 || n1>=s.Nu){
			return Ddata[p][n2][q2];
		}else if(n2<0 || n2>=s.Nu){
			return Ddata[p][n1][q2];
		}else{
			return wu1*Ddata[p][n1][q2] + wu2*Ddata[p][n2][q2];
		}
		
	}else if(q2<0 || q2>=s.Nv){
		if(n1<0 || n1>=s.Nu){
			return Ddata[p][n2][q1];
		}else if(n2<0 || n2>=s.Nu){
			return Ddata[p][n1][q1];
		}else{
			return wu1*Ddata[p][n1][q1] + wu2*Ddata[p][n2][q1];
		}
		
	}else{
		if(n1<0 || n1>=s.Nu){
			return wv1*Ddata[p][n2][q1] + wv2* Ddata[p][n2][q2];
		}else if(n2<0 || n2>=s.Nu){
			return wv1*Ddata[p][n1][q1] + wv2*Ddata[p][n1][q2];
		}else{
			return wv1*(wu1*Ddata[p][n1][q1] + wu2*Ddata[p][n2][q1] ) 
		+ wv2*(wu1*Ddata[p][n1][q2] + wu2*Ddata[p][n2][q2] );
		}	
	}
}

void backproj(ScanData & s)
// sets the backprojected data.
{
	int p=0; // t index
	int n=0; // u index
	int q=0; // v index
	int i=0; // x index
	int j=0; // y index
	int k=0; // z index
	int pi,pf,p1;
	Float x=0,y=0,z=0,G=0,ta=0,tb=0, u=0, v=0;

	vector< vector< vector<Float> > >  Ddata (s.Nt, vector< vector<Float> > (s.Nu, vector<Float> (s.Nv,0)));

	// Differentiate
	//u derivative
	for(p=0; p<s.Nt;++p){ // t loop
		for(q=0; q<s.Nv; ++q){ // v loop
			v = s.vmin + s.dv*q;
			u = s.umin;
			Ddata[p][0][q] = (s.D*s.D + u*u)/s.D*(-1.5f*s.data[p][0][q] + 2.0f*s.data[p][1][q] - 0.5f*s.data[p][2][q])/s.du;
			for(n=1; n<s.Nu-1;++n){ // u loop
				u = s.umin + s.du*n;
				Ddata[p][n][q] = (s.D*s.D + u*u)/s.D*(-0.5f*s.data[p][n-1][q] + 0.5f*s.data[p][n+1][q])/s.du;
			}
			u = s.umax;
			Ddata[p][s.Nu-1][q] = (s.D*s.D + u*u)/s.D*(1.5f*s.data[p][s.Nu-1][q] - 2.0f*s.data[p][s.Nu-2][q] + 0.5f*s.data[p][s.Nu-3][q])/s.du;
		}
	}

	//v derivative
	for(p=0; p<s.Nt;++p){ // t loop
		for(n=0; n<s.Nu; ++n){ // u loop
			u = s.umin + s.du*n;
			v = s.vmin;
			Ddata[p][n][0] += u*v/s.D*(-1.5f*s.data[p][n][0] + 2.0f*s.data[p][n][1] - 0.5f*s.data[p][n][2])/s.dv;
			for(q=1; q<s.Nv-1;++q){ // v loop
				v = s.vmin + s.dv*q;
				Ddata[p][n][q] += u*v/s.D*(-0.5f*s.data[p][n][q-1] + 0.5f*s.data[p][n][q+1])/s.dv;
			}
			v = s.vmax;
			Ddata[p][n][s.Nv-1] += u*v/s.D*(1.5f*s.data[p][n][s.Nv-1] - 2.0f*s.data[p][n][s.Nv-2] + 0.5f*s.data[p][n][s.Nv-3])/s.dv;
		}
	}
	
	//t derivative
	for(q=0; q<s.Nv; ++q){
		for(n=0; n<s.Nu;++n){ // u loop
			Ddata[0][n][q] = (-1.5f*s.data[0][n][q] + 2.*s.data[1][n][q] - 0.5f*s.data[2][n][q])/s.dt - Ddata[0][n][q];
			for(p=1; p<s.Nt-1;++p){ // t loop
				Ddata[p][n][q] = (-0.5f*s.data[p-1][n][q] + 0.5f*s.data[p+1][n][q])/s.dt - Ddata[p][n][q];
			}
			Ddata[s.Nt-1][n][q] = (1.5f*s.data[s.Nt-1][n][q] - 2.0f*s.data[s.Nt-2][n][q] + 0.5f*s.data[s.Nt-3][n][q])/s.dt - Ddata[s.Nt-1][n][q];
		}
	}

	// Preweight
	for(p=0;p<s.Nt;++p){ // t loop
		for(n=0; n<s.Nu; ++n){ // u loop
			u = s.umin + s.du*n;
			for(q=0; q<s.Nv; ++q){ // v loop
				v = s.vmin + s.dv*q;
				Ddata[p][n][q] = Ddata[p][n][q] * s.D / sqrt(s.D*s.D + u*u + v*v);
			}
		}
	}

	// Perform backprojection
	for(j=0; j<s.Ny;++j){ // y loop
		y = s.ymin + s.dy * j;
		ta = ShiftPosition( asin(y/s.R) );
		tb = ShiftPosition( PI - ta );
		pi = (int) PositionIndex(s ,ta,true);
		pf = (int) PositionIndex(s ,tb,false);

		for(k=0; k<s.Nz;++k){
		//for(k=255;k<=256;k++){ // central slice only
			z = s.zmin + s.dz*k;
			for(i=0; i<s.Nx;++i){ // x loop
				x = s.xmin + s.dx * i;
				if(s.DBPmask[i][j][k]){/* to be uncommented when project supports conjugate rays
					if(isMissing(s,ta)){
						project(s, Ddata, s.DBPdata[i][j], 0,  pf,  x,  y, z, false);
						pi = (int) ceil((ta-s.tmin)/s.dt);
						project(s, Ddata, s.DBPdata[i][j], pi,  -1,  x,  y, z, true);
					}else if(isMissing(s,tb)){
						project(s, Ddata, s.DBPdata[i][j], pi,  s.Nt-1,  x,  y, z, false);
						if(tb<0) tb = tb + 2.* PI;
						pf = (int) ceil((tb-s.tmax)/s.dt);
						project(s, Ddata,s.DBPdata[i][j], s.Nt,  pf+s.Nt-1,  x,  y, z, true);
					}else{
						if(pi>pf){
							project(s, Ddata, s.DBPdata[i][j], pi,  s.Nt-1,  x,  y, z, false);
							p1 = (int) floor((s.tmin-s.tmax)/s.dt); //might need to be improved
							project(s, Ddata, s.DBPdata[i][j], s.Nt,  s.Nt+p1-1,  x,  y, z, true);
							project(s, Ddata, s.DBPdata[i][j], 1,  pf,  x,  y, z, false);
						
						}else{*/
							project(s, Ddata, s.DBPdata[i][j][k], ta, tb, pi,  pf,  x,  y, z, false);
						//}
					//}
				}
			}
		}
	}
}

Float atan3(Float y,Float x)
{
	Float a = atan2(y,x);
	if (a>=0){ return a;}
	else{ return 2.*PI+a;}
}

Float ShiftPosition(const Float t){
	Float T = fmod((double) t + 9.*PI,2.*PI);
	return T - PI;
}

Float DetectorPosition(Float x, Float y, Float t, Float R){
	Float phi = atan3(y-R*sin(t) , x-R*cos(t) );
	if( t<=PI && t>=PI/2 && phi < PI ){
		return phi - t + PI;
	}else if(t>-PI && t<=-PI/2 && phi>=PI ){
		return phi - t - 3.* PI;
	}else{
		return phi - t - PI;
	}
}

Float IndexPosition(const ScanData &s, const int p){
	return ShiftPosition(s.tmin + s.dt * p);
}

int PositionIndex(const ScanData &s, const Float t, bool isCeil){
	Float T =(t-s.tmin);
	if ((ShiftPosition(t)<s.tmin||ShiftPosition(t)>s.tmin+PI)&&!(s.tmin<0 && t>ShiftPosition(PI+s.tmin))) T= 2.*PI+ T;
	
	int p;
	if(isCeil){p = (int) ceil(fmod((double)(T)/s.dt + 8.*(s.Nt),s.Nt));}
	else{p = (int) floor(fmod((double)(T)/s.dt + 8.*(s.Nt),s.Nt));}
	return p;
}

bool isMissing(const ScanData & s, const Float t)
{
	if(s.tmin < 0 && s.tmax < 0){
		if(s.tmin > s.tmax && t>s.tmax && t<s.tmin){ return true;}
		else if(s.tmin < s.tmax && (t>s.tmax || t<s.tmin)){ return true;}
		else {return false;}
	}else if(s.tmin >= 0 && s.tmax >= 0){
		if(s.tmin > s.tmax && t>s.tmax && t<s.tmin){ return true;}
		else if(s.tmin < s.tmax && (t>s.tmax || t<s.tmin)){ return true;}
		else {return false;}
	}else if(s.tmin>=0 && s.tmax<0){
		if(t>s.tmax && t<s.tmin){return true;}
		else{return false;}
	}else if(s.tmin< 0 && s.tmax>=0){
		if(t>s.tmax || t<s.tmin){return true;}
		else{return false;}
	}else{
		return false;
	}
}

bool isHilbertKnown(const ScanData & s, Float x, Float y, Float z)
{
	if (x*x + y*y <= s.rov*s.rov){ //return true;
		Float medang = (s.tmin + s.tmax)/2.;
		if(s.tmax>s.tmin) medang = medang + PI;

		Float k[2] = {cos(medang), sin(medang)};
		Float C = k[0]*s.R*s.cost[s.Nt-1] + k[1]*s.R*s.sint[s.Nt-1];

		if(k[0]*x + k[1]*y < C){
			Float zcut = (s.R - sqrt(x*x + y*y))*s.vmax/s.D;
			if (abs(z) < zcut){ return true;}
			else{return false;}
		}else{return false;}
	}else{
		return false;
	}
}

void project(ScanData & s, const vector< vector< vector<Float> > > & Ddata, Float & image, Float ta, Float tb, int pi, int pf, Float x, Float y, Float z, bool isConj){
	int p=0,p1=0;
	Float u, v, L, U, value, ti, tf, di, df,w_i0, w_i1, w_f1, w_f0;
	ti = IndexPosition(s,pi);
	tf = IndexPosition(s,pf);
	di = ti-ta;
	df = tb-tf;
	if(di<0 || di>s.dt ||df<0 || df>s.dt ) cout << di << " " <<df<<endl;
	w_i1 = di*di*0.5/s.dt;
	w_i0 = 0.5*di*(2.0-di/s.dt) + 0.5*s.dt;
	w_f1 = df*df*0.5/s.dt;
	w_f0 = 0.5*df*(2.0-df/s.dt) + 0.5*s.dt;


	if(isConj){
		for(p=pi;p<=pf;p++){/* To be modified to get conjugate ray
			t = ShiftPosition(s.tmin + p*s.dt);
			L = sqrt( pow(x - s.R*s.cost[p],2) + pow(y - s.R*s.sint[p],2) );
			G = DetectorPosition(x, y, t, s.R);			
			t = ShiftPosition(t + 2.*G + PI);
			p1 = PositionIndex(s ,t,true);
			image = image - s.dt/L * DFint(s, Ddata, p1,-G);*/
		}
	}else{
		for(p=pi-1; p<=pf+1;p++){
			//L = sqrt(pow(x-s.R*s.cost[p] , 2) + pow(y-s.R*s.sint[p],2));
			U = s.R - x*s.cost[p] - y*s.sint[p]; // To be checked
			u = s.D/U*(x*s.sint[p]-y*s.cost[p]);
			v = s.D/U*z;
			//U=L;
			if (p == pi){image += (w_i0)/U *DFint(s, Ddata, p, u, v);}
			else if (p == pi-1){image += (w_i1)/U * DFint(s, Ddata, p, u, v);}
			else if(p == pf){image += (w_f0)/U * DFint(s, Ddata, p, u, v);}
			else if(p == pf+1){image += (w_f1)/U * DFint(s, Ddata, p, u, v);}
			//else if(p == pi-1){image += (w_i1)/U * DFint(s, Ddata, p, u, v);}
			//else if(p == pi+1){image += (w_f1)/U * DFint(s, Ddata, p, u, v);}
			else{image += s.dt/U * DFint(s, Ddata, p, u, v);}
			//image += s.dt/U * DFint(s, Ddata, p, u, v);
			//image = image + s.dt * s.R*s.D/(U*U) * DFint(s, Ddata, p, u, v);
		}
	}
}



void backproj_ibp(ScanData & s)
// sets the backprojected data.
{
	int p=0; // t index
	int n=0; // u index
	int q=0; // v index
	int i=0; // x index
	int j=0; // y index
	int k=0; // z index
	int pi,pf,p1;
	Float x=0,y=0,z=0,G=0,ta=0,tb=0, u=0, v=0,U,bterma1,bterma2,btermb1,btermb2, ta1,ta2,tb1,tb2;

	vector< vector< vector<Float> > >  Ddata (s.Nt, vector< vector<Float> > (s.Nu, vector<Float> (s.Nv,0)));
	vector< vector< vector<Float> > >  Wdata (s.Nt, vector< vector<Float> > (s.Nu, vector<Float> (s.Nv,0)));


	// Preweight
	for(p=0;p<s.Nt;++p){ // t loop
		for(n=0; n<s.Nu; ++n){ // u loop
			u = s.umin + s.du*n;
			for(q=0; q<s.Nv; ++q){ // v loop
				v = s.vmin + s.dv*q;
				Wdata[p][n][q] = s.data[p][n][q] * s.R / sqrt(s.D*s.D + u*u + v*v);
			}
		}
	}

	// Differentiate
	//u derivative
	for(p=0; p<s.Nt;++p){ // t loop
		for(q=0; q<s.Nv; ++q){ // v loop
			v = s.vmin + s.dv*q;
			u = s.umin;
			Ddata[p][0][q] = (-1.5f*Wdata[p][0][q] + 2.0f*Wdata[p][1][q] - 0.5f*Wdata[p][2][q])/s.du;
			for(n=1; n<s.Nu-1;++n){ // u loop
				u = s.umin + s.du*n;
				Ddata[p][n][q] = (-0.5f*Wdata[p][n-1][q] + 0.5f*Wdata[p][n+1][q])/s.du;
			}
			u = s.umax;
			Ddata[p][s.Nu-1][q] = (1.5f*Wdata[p][s.Nu-1][q] - 2.0f*Wdata[p][s.Nu-2][q] + 0.5f*Wdata[p][s.Nu-3][q])/s.du;
		}
	}


	// Perform backprojection
	for(j=0; j<s.Ny;++j){ // y loop
		y = s.ymin + s.dy * j;
		ta = ShiftPosition( asin(y/s.R) );
		tb = ShiftPosition( PI - ta );
		pi = (int) PositionIndex(s ,ta,true);
		pf = (int) PositionIndex(s ,tb,false);
		
		for(k=0; k<s.Nz;++k){
		//for(k=255;k<=256;k++){ //for central slice only
			z = s.zmin + s.dz*k;
			for(i=0; i<s.Nx;++i){ // x loop
				x = s.xmin + s.dx * i;
				if(s.DBPmask[i][j][k]){
					U = s.R - x*s.cost[pi] - y*s.sint[pi]; 
					u = s.D/U*(x*s.sint[pi]-y*s.cost[pi]);
					v = s.D/U*z;
					bterma1 = DFint(s, s.data, pi, u, v)/sqrt(pow(x - s.R*s.cost[pi],2) + pow(y - s.R*s.sint[pi],2));

					U = s.R - x*s.cost[pi-1] - y*s.sint[pi-1];
					u = s.D/U*(x*s.sint[pi-1]-y*s.cost[pi-1]);
					v = s.D/U*z;
					bterma2 = DFint(s, s.data, pi-1, u, v)/sqrt(pow(x - s.R*s.cost[pi-1],2) + pow(y - s.R*s.sint[pi-1],2));

					U = s.R - x*s.cost[pf] - y*s.sint[pf]; 
					u = s.D/U*(x*s.sint[pf]-y*s.cost[pf]);
					v = s.D/U*z;
					btermb1 = DFint(s, s.data, pf, u, v)/sqrt(pow(x - s.R*s.cost[pf],2) + pow(y - s.R*s.sint[pf],2));

					U = s.R - x*s.cost[pf+1] - y*s.sint[pf+1]; 
					u = s.D/U*(x*s.sint[pf+1]-y*s.cost[pf+1]);
					v = s.D/U*z;
					btermb2 = DFint(s, s.data, pf+1, u, v)/sqrt(pow(x - s.R*s.cost[pf+1],2) + pow(y - s.R*s.sint[pf+1],2)); 

					ta1 = IndexPosition(s,pi);
					ta2 = IndexPosition(s,pi-1);
					tb1 = IndexPosition(s,pf);
					tb2 = IndexPosition(s,pf+1);

					s.DBPdata[i][j][k] = bterma1*(ta-ta2)/(ta1 - ta2) + bterma2*(ta-ta1)/(ta2 - ta1) - btermb1*(tb-tb2)/(tb1 - tb2) - btermb2*(tb-tb1)/(tb2 - tb1);
					project_ibp(s, Ddata, s.DBPdata[i][j][k], ta, tb, pi,  pf,  x,  y, z, false);

				}
			}
		}
	}
}


void project_ibp(ScanData & s, const vector< vector< vector<Float> > > & Ddata, Float & image, Float ta, Float tb, int pi, int pf, Float x, Float y, Float z, bool isConj){
	int p=0,p1=0;
	Float u, v, L, U, value, ti, tf, di, df,w_i0, w_i1, w_f1, w_f0;
	ti = IndexPosition(s,pi);
	tf = IndexPosition(s,pf);
	di = ti-ta;
	df = tb-tf;
	if(di<0 || di>s.dt ||df<0 || df>s.dt ) cout << di << " " <<df<<endl;
	w_i1 = di*di*0.5/s.dt;
	w_i0 = 0.5*di*(2.0-di/s.dt) + 0.5*s.dt;
	w_f1 = df*df*0.5/s.dt;
	w_f0 = 0.5*df*(2.0-df/s.dt) + 0.5*s.dt;


	if(isConj){
		for(p=pi;p<=pf;p++){/* To be modified to get conjugate ray
			t = ShiftPosition(s.tmin + p*s.dt);
			L = sqrt( pow(x - s.R*s.cost[p],2) + pow(y - s.R*s.sint[p],2) );
			G = DetectorPosition(x, y, t, s.R);			
			t = ShiftPosition(t + 2.*G + PI);
			p1 = PositionIndex(s ,t,true);
			image = image - s.dt/L * DFint(s, Ddata, p1,-G);*/
		}
	}else{
		for(p=pi-1; p<=pf+1;p++){
			//L = sqrt(pow(x-s.R*s.cost[p] , 2) + pow(y-s.R*s.sint[p],2));
			U = s.R - x*s.cost[p] - y*s.sint[p]; // To be checked
			u = s.D/U*(x*s.sint[p]-y*s.cost[p]);
			v = s.D/U*z;
			//U=L;
			if (p == pi){image += (w_i0)*s.D*s.D/U/U *DFint(s, Ddata, p, u, v);}
			else if (p == pi-1){image += (w_i1)*s.D*s.D/U/U * DFint(s, Ddata, p, u, v);}
			else if(p == pf){image += (w_f0)*s.D*s.D/U/U * DFint(s, Ddata, p, u, v);}
			else if(p == pf+1){image += (w_f1)*s.D*s.D/U/U * DFint(s, Ddata, p, u, v);}
			//else if(p == pi-1){image += (w_i1)/U * DFint(s, Ddata, p, u, v);}
			//else if(p == pi+1){image += (w_f1)/U * DFint(s, Ddata, p, u, v);}
			else{image += s.dt*s.D*s.D/U/U * DFint(s, Ddata, p, u, v);}
			//image += s.dt/U * DFint(s, Ddata, p, u, v);
			//image = image + s.dt * s.R*s.D/(U*U) * DFint(s, Ddata, p, u, v);
		}
	}
}
